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Abstract 

We explore the rotational degree of freedom between graphene layers via the simple prototype 
of the graphene twist bilayer, i.e., two layers rotated by some angle 6. It is shown that, due to 
the weak interaction between graphene layers, many features of this system can be understood by 
interference conditions between the quantum states of the two layers, mathematically expressed as 
Diophantine problems. Based on this general analysis we demonstrate that while the Dirac cones 
from each layer are always effectively degenerate, the Fermi velocity vp of the Dirac cones decreases 
as 6 — > 0°; the form we derive for vf{9) agrees with that found via a continuum approximation in 
Ref. From tight binding calculations for structures with 1.47° < 6 < 30° we find agreement 
with this formula for 9 > 5°. In contrast, for 9 < 5° this formula breaks down and the Dirac bands 
become strongly warped as the limit 6 — > is approached. For an ideal system of twisted layers 
the limit as 9 — > 0° is singular as for 9 > the Dirac point is fourfold degenerate, while at 9 = 
one has the twofold degeneracy of the AB stacked bilayer. Interestingly, in this limit the electronic 
properties are in an essential way determined globally, in contrast to the 'nearsightedness— of 
electronic structure generally found in condensed matter. 

PACS numbers: 73.20.At, 73.21.Ac, 81.05.Uw 
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I. INTRODUCTION 



In addition to offering a possible route towards exploiting the many remarkable properties 
of graphene^, the epitaxial growth of graphene on SiC^~- presents a number of mysterious 
aspects. Principle amongst these is that the thermally induced growth of graphene on 
the C-face typically results in several graphene layers and yet, remarkably, this complex 
graphene-based system shows behavior identical to that of single layer graphene (SLG). In 
striking contrast, bilayer graphene produced by mechanical exfoliation has already a different 
low energy electronic structure to that of SLG; a quadratic dispersion instead of linear. 

An insight into this intriguing behavior of the C-face growth was recently provided by 
Hass et al.—. These authors showed that growth on the C face results in a high density of 
twist boundary faults, i.e., layers with a relative rotation. Furthermore, ab-initio calcula- 
tions by the same authors showed that if two graphene layers were rotated with the same 
relative rotation observed in experiment, 9 = 30° ± 2.20°, then these layers exhibited a lin- 
ear spectrum near the Dirac point, exactly as in SLG. Rotation and translation of graphene 
layers thus have profoundly different impact on the low energy spectrum, and this lies at 
the heart of the C-face behavior. 

While rationalizing the SLG nature of the C-face, these findings raised a number of 
questions. Firstly, as to the character of the rotational degree of freedom in few layer 
graphene systems: do all rotations cause such an electronic decoupling or, alternatively, 
only a subset of "magic" angles? This question is relevant to experiments as subsequent 
investigations have shown that various angles of rotation may occur during growth on the 
C-face^. Clearly, a related question is the nature of the mechanism responsible for this 
electronic decoupling: how does the rotation lead to the emergence of an effective Dirac- 
Weyl equation for low energies? 

These questions, at first sight, appear difficult from the point of view of theory as one 
ultimately requires general statements to be made about an infinite class of possible lattices. 
Initially, theoretical progress was made by example of specific rotation angles or limits, with 
graphene bilayer and trilayer systems calculated ab-initio in Ref. 10(, while in Ref. |l| the 
— > 0° limit of the twist bilayer was investigated via a continuum approximation to the 
tight binding Hamiltonian. In the former case a low energy linear spectrum was noted for 
all layers experiencing a relative rotation, while the latter work found also a linear spectrum 
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but with the Fermi velocity at the Dirac point, Vf, strongly suppressed as compared to SLG. 
Subsequent Raman spectroscopy experiments^^ differ on whether this effect is present in 
misoriented graphene samples; in Ref. ll| a blueshift of the graphene 2D peak was attributed 
to this effect, however in Ref. [12j this was instead attributed to a modification of the phonon 
dispersion in misoriented layers. 

In Ref. 13|] it was shown that the rotational degree of freedom was associated with a de- 



structive interference of quantum states from each layer, and that this resulted in a coupling 
that becomes progressively weaker as the size of the commensuration cell increases. In fact, 
coupling at the Dirac point is already very weak for the smallest possible commensuration, 
a cell of 28 carbon atoms, with a splitting of 7meV found in ab-initio calculations^. All 
misoriented graphene layers are, therefore, predicted to show effectively decoupled Dirac 
cones. 

Further theoretical investigations have been undertaken with regard to both the energetics 
of misoriented layers^, and the simulation of scanning tunneling microscopy images for 
such layers^. In the former work it was noted that the sliding energy of relatively rotated 
graphene layers is essentially zero, in dramatic contrast to the case without rotation where 
the AB configuration is energetically favored. Most recently, tight-binding calculations have 
been performed for a wide range of misorientation angles^. This latter work demonstrates 
a reduction in the Fermi velocity that, for a wide range of rotation angles, agrees with the 
result of Ref. fl. 

In this article we aim to accomplish two things. Firstly, the formalism presented in 
Ref. [fij is extended to explain, on general lattice grounds (i.e., without deploying a con- 
tinuum approximation), both the Dirac cone decoupling and Fermi velocity suppression. 
Secondly, we provide a numerical implementation of this formalism using the tight-binding 
method. We demonstrate that this numerical scheme is at least an order of magnitude 
faster than the usual tight-binding basis, and using this explore the electronic structure as 
a function of rotation angle for 1.47° < 9 < 30°. 

We now present a brief summary of the content of this article. Firstly, in Section II 
we discuss in detail the crystal structure of mutually rotated graphene layers, and derive 
the conditions for a commensurate crystal structure to occur. An important feature of this 
system is the emergence, for 9 < 15°, of a so-called moire pattern^. This is a hexagonal 
interference pattern, consisting of regions of AA and AB stacking, the periodicity of which 
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represents a new structural length scale of the system. 

Section III then describes the electronic structure of the bilayer in terms of a basis formed 
from the quantum states of the two mutually rotated layers with, in addition, the bilayer one- 
electron potential treated as a superposition of two single layer potentials, i.e. + R^ 2 ** 
with R the rotation operator and V^ 1 ' 2 ' potentials with the in-plane translation symmetry 
of SLG, an approach first described in Ref. [lsj]. It is shown how this leads to a convenient 
separating out of purely symmetry related aspects of the electronic structure, leading to 
simple conditions for determining if the overlap elements of the potential with single layer 
states are vanishing or not. As the interlayer interaction part of the full bilayer Hamiltonian 
may be constructed from such overlap elements, an understanding of how and why these 
vanish leads in turn to an understanding of the nature of the interlayer decoupling in this 
system. 

In this context we investigate how the overlap between states from the constituent layers 
depends on their k- vectors (i.e., their k- vectors in the two mutually rotated single layer 
Brillouin zones). We find that this dependence is rather subtle, and that the vanishing or 
not of such overlaps depends crucially on these k-vectors. On this basis we demonstrate 
a number of general features of the bilayer electronic structure, and in particular show 
that (i) for the Dirac bands the 1st order term of a perturbation theory in the interlayer 
interaction is negligible for all rotations and that, furthermore, (ii) 2nd order order terms in 
perturbation theory lead to a Fermi velocity suppression of the form found in Ref. We 
further develop two consequences of (i); if 2nd (and higher) order terms are unimportant 
then the Dirac cones effectively decouple, and that the Dirac bands from each layer are 
degenerate regardless of the role of higher order terms. 

Discussed also in this Section is the rather unusual 9 — > limit, which is a singular 
limit as for any 9 > the electronic structure is dramatically different from that at 9 = 0. 
This is, of course, simply a electronic manifestation of the fact that the lattice geometry is 
also singular in this limit: for any small but non-zero 9 one has a moire pattern, while at 
9 = the graphene layers are simply AB (or AA) stacked. An interesting aspect of this 
limit is a breakdown in the notion of 'nearsightedness', i.e., that electronic properties are 
essentially determined locally. As the moire periodicity diverges as 9 — > 0° and, furthermore, 
as the twist bilayer electronic structure must, for any finite 9, be different from both the 
AA and AB stacked bilayers, one concludes that in this limit the electronic properties are, 
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in contrast, in an essential way determined globally. 

Finally, Section IV is devoted to a presentation of tight-binding calculations for the 
graphene twist bilayer. We demonstrate that a basis formed by the quantum states of the two 
mutually rotated layers converges remarkably quickly, and leads to a dramatic improvement 
in computational efficiency Using this we then investigate the bilayer electronic structure 
for 1.47° < 9 < 30° and find a suppression of the Fermi velocity, vf, that is dramatic for 
small angles (at 9 = 1.47° the reduction in vf is 95%) but, in agreement with all ab-initio 
calculations to date^^ 1 ^, insignificant for 9 >m 15°. However, while the expression for 

n 

the Fermi velocity suppression derived here and in Ref. [1], describes almost perfectly the 
tight-binding results for 5° < 9 < 30°, it breaks down for 9 < 5°. This breakdown is a result 
of the fact that the Fermi velocity suppression is a 2nd order effect in layer interaction, and 
in the limit 9 — > 0° the resulting near degeneracy of the Dirac cones entails the importance 
of terms beyond this order. 



II. COMMENSURATION CONDITIONS OF THE TWISTED BILAYER 



A prerequisite to exploring the electronic structure of the twisted bilayer is an elucidation 
of the crystallography of such a system, i.e., determining the conditions under which two 
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misoriented layers are in commensuration. This problem was studied in Ref. 
complete solution was presented; here we provide a more detailed derivation of those results 
with, additionally, a somewhat simpler and more symmetrical choice of commensuration 
vectors. 

Evidently, the existence of a commensuration depends only on the relative rotation of the 
lattice vectors of each layer, and not on the structure of the unit cells of each layer. Thus 
we need not, at this stage, concern ourselves with which axis the rotation is taken about 
and the initial configuration (AB or AA, and so on) of the graphene bilayer; these amount 
to different choices of initial basis vectors within each cell. The commensuration condition 
may be written as 



ri = Rr 2 (1) 

where ri, r 2 are hexagonal lattice vectors, and R the rotation operator. The set {ri} are the 
resulting coincident points between the two layers, while {r 2 } is the same set, but viewed 
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FIG. 1: (Color online) Shown is the number of C atoms, Nc, in the commensuration cell as a 
function of the relative orientation of the two graphene layers, for Nc < 7000. Inset displays the 
moire pattern for the cell indicated number 4. Band structures of twist bilayers corresponding to 
the points labeled 1-4 are displayed in panels 1-4 of Fig. \7\ The dashed line corresponds the lower 
bound Nc = (sin 2 0/2)~ 1 ; for commensuration cells that fall on this line the moire periodicity is 
equal to the commensuration periodicity, see Section II for details. 
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from the rotated coordinate system. Utilizing the unrotated lattice as a coordinate system, 
i.e., r = ia 1 +ja 2 with i,j E Z, Eq. (JTJ, with a standard choice of lattice vectors a x = [V3, 0] 
and &2 = [^, 3/2], may be written as 



/ mi ' 




\m 2 J 





cos# — ^sin^ — ^sin^ 



cos ^ + ^75 sm @ 




(2) 



(Here and throughout this article our unit of length is chosen to be the graphene C-C 
separation.) This maps one integer pair (711,712) to another (7771,7712) and, for this to be 
possible, a necessary and sufficient condition on the matrix in Eq. (J2]) is that it assumes only 
rational values^. This leads to the following conditions on 9 
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FIG. 2: (Color online) Illustration of the commensuration cell for the case of a misorientation angle 
of 8 = 21.78°, generated by a (p, q) pair of (1,3) (lattice vectors ti and t2). Shown also are the unit 
cells of the unrotated graphene layer (vectors ai, a2) and rotated graphene layer (vectors Rai, 
Ra2). For explanations of other symbols refer to Section 2. 




-^sinfl = l ± (3) 

cos9 = ^, (4) 

where ii, «2) *3 £ N. These are therefore related by the following second order Diophantine 
equation 



3il + i 2 2 = il (5) 

Solution of this equation proceeds in a standard way-^ (analogous to the case of Pythagorean 
triples) by dividing by i\ and making the substitution x — ^, y = ^. There is thus a one 
to one mapping between solutions of Eq. (JHJ) and rational points on the ellipse 3x 2 + y 2 = 1. 
One such point is (0, 1) and any other may be found by the intersection with the ellipse of 
a line passing through (0, 1) and (q/p,0). The coordinates of this latter point then lead to 
the following solution for ii, i 2 , «3 
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z'l = 2pq (6) 
12 = 3g 2 - p 2 (7) 
*s = 3g 2 +p 2 (8) 

where p, q G N. From these equations we immediately find the set of rotation angles leading 
to commensurations, 

For q > p > this formula produces rotation angles that lie in the range < < 60°. All 
other rotation angles are equivalent due to the symmetry of the hexagonal lattice. Clearly, 
the limit p/q — > corresponds to — > 0° while, on the other hand, the limit p/q — > 1 
corresponds to — > 60°. Note that changing the sign of p or q sends — > —0. Since the 
limit — > 60° is equivalent to — > 0° taken from below, and since formally little changes by 
the substitution p — > — p or g — >■ — g, our focus in this work will be on angles in the range 
0° < < 30°. 

We also require the corresponding primitive vectors of the commensuration lattice. Sub- 
stitution of Eqs. 03]) and (jl]) into Eq. (T5]) results in the following coupled linear Diophantine 
equations 



/ rrii ' 




\m 2 J 






(10) 



i 2 - i\ -2ii 
2ii i 2 + i\ 

The solution of these equations follows by a similarity transform such that the matrix mul- 
tiplying {n\n,2) T is diagonal. Crucially, the eigenvectors of this matrix are independent of 
p, q and thus the problem is recast as coupled linear diophantine equations that are linear 
in p, q. These may then be solved by inspection yielding the result that 
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TABLE I: Possible values that the parameter 7 can take. 





5 = 1 


(5 = 3 


p, q odd 


6 


2 


otherwise 


3 


1 



with a, (3 are arbitrary constants such that n and m are integer valued. The final step is to 
determine the primitive vectors of the commensuration lattice. This calculation we present 
in Appendix A, and here quote only the result. The form of the commensuration vectors 
turns out to depend on a parameter 5 = 3/gcd(p, 3). For the case where 5 = 1 we find 



1 / P + 3<? 
7 \ -2p 




(13) 



while for the case 5 = 3 we find 



*, = i ~ p - q U=M 29 ], (14) 

-y\ 2q ) 1 \-p + q' 

where 7 = gcd(3g + p, 3q — p). Values that this parameter may take when gcd(p, q) — 1 are 
indicated in Table [H 

Thus every possible commensuration between misoriented layers is uniquely specified by 
an integer pair p > q > such that gcd(p, q) = 1. Given this we can completely characterize 
the commensuration; the rotation angle may be obtained from Eq. while the lattice 
vectors are given by either Eq. ffTB"]) or Eq. (JHJ) , depending on whether the parameter 
6 = 3/gcd(p, 3) assumes the values of 1 or 3 respectively. The various notations introduced 
in this derivation are illustrated in Fig. [2j 

It is worth reflecting on the reason that two integers, p and q, are needed to specify a 
commensuration while, on the other hand it is clear that any bilayer lattice (commensurate 
or incommensurate) is uniquely specified by a single number 9. This is a consequence of 
the relation between the real and rational number fields: given a 9 there are infinitely many 
choices of p and q in Eq. [9] such that 9 may be reproduced to an arbitrary accuracy e. 

The ratio of the total number of lattice vectors to coincident lattice vectors for the twist 
boundary is found to be given by 
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with the number of carbon atoms in the commensuration cell Nc = 4JV. (The factor 4 
simply arising from the fact there are two layers in the cell, and two basis atoms in the 
honeycomb structure.) In Fig. [TJis plotted Nc, as a function of misorientation angle; the 
minimum Nq is 28 corresponding to 9 = 30° ± 8.21, however Nc diverges in the 9 — > 0° (or 
9 — > 60°) limits. Combining Eqs. © and (fT5|) we find that 

N= Ts^m- (16) 

and so for 9 0° N diverges as 1/9 2 . This small angle limit is associated with the emergence 
of a new structural length scale, that of the moire periodicity D; such a moire pattern is 
shown in the inset of Fig. [TJ The relation between D and 9 is given by^ 

where a is the graphene lattice constant. The relation between the lattice constant of the 
commensuration cell and the moire periodicity may be seen by setting p = 1, 8 = 3, 7 = 2 
in Eq. (116I) . (corresponding to cells generated by p = 1 and q an odd integer, see Tabled]), 
and using N = D 2 /a 2 which then gives back the formula for the moire periodicity, Eq. ( TT7I) . 
In this case, therefore, the moire periodicity is equal to the commensuration cell lattice 
constant. For these commensuration cells we find Nc = (sin 2 #/2) -1 , and this is the lower 
bound function plotted in Fig. [lj On the other hand, for all other commensuration cells the 
"commensuration periodicity" is greater than the moire periodicity. 

Finally, we note that the analytic results presented here are in agreement with the nu- 
merical solution to this problem provided recently by Campenara et a/.—; special cases of 
these results have been found in Ref. [l| (the case p = 1, 5 = 3, 7 = 2) and more recently in 



Ref. 



16| (the case S = 1). 



III. ANALYSIS OF THE INTERLAYER INTERACTION 



In this section we shall describe how the problem of understanding the general electronic 



properties of the twist bilayer for any 9 is solved. Our approach is that described in Ref. 
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FIG. 3: (Color online) Brillouin zones of the unrotated (U) and rotated (R) graphene layers, as 
well as the Brillouin zone of the bilayer supercell (B) for the case (p,q) = (1,5), corresponding 
to 9 = 13.17°. In this Figure b>i, l>2 are the reciprocal lattice vectors of the unrotated graphene 
layer, Rbi, Rb2 the reciprocal lattice vectors of the rotated layer, and gi, g2 the reciprocal lattice 
vectors of the bilayer supercell. Special K points of various Brillouin zones indicated by subscript 
U, R, and B. The separations of special K points indicated are AK = |Kr — K[/| = |K|j — K^|, 
Aifi = \K*rj - K R \, and AK 2 = \K V - KU. 
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which, in broad outline, may be characterized as 'constructing the bilayer system from 
single layer components'. We take the bilayer potential as a superposition of single layer 
potentials, and use as a basis for this new system the eigenkets of the single layer systems. 
The advantage of this is that the resulting matrix elements may then be analyzed as a 
commensuration problem of reciprocal space lattices. Such commensuration problems can 
be readily understood for any angle, and thus one may then understand the physics of the 
twist bilayer for general angle. 

The remainder of this Section proceeds as follows. In Section III A we first set up the 
Hamiltonian and basis used to analyze the twist bilayer. Following this, in Section III B the 
various reciprocal lattices and associated Brillouin zones are described. In Section III C a 
condition is derived that determines the vanishing of overlap elements found in the model of 
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Section III A. Finally, in Sections III D-F, we use this understanding to determine a number 
of generic electronic properties of the twist bilayer. 



A. Model Hamiltonian and basis 

In considering the interaction between misoriented layers it is useful to take the bilayer 
potential as simply the superposition of two single layer potentials, i.e., 

H=^ + V«+V( 2 ). (18) 
2m 

Here V"( n ) are one-electron single layer graphene (SLG) potential operators that satisfy 



H (n) 



2 



(n)\ _ » 



(19) 



where are SLG eigenvalues, and % and k represent band and k-vector quantum num- 
bers respectively. These one-electron SLG potentials are invariant under different in-plane 
translations; we have TV* 1 ) = V* 1 ) and RTR _1 V (2 ^ = V* 2 ). We shall take the superscript 
"1" to denote objects associated with the unrotated layer, and superscript "2" for objects 
associated with the rotated layer. (Such a designation is, in mutually rotated layers, clearly 
arbitrary and made only for convenience.) Given the weak interaction (and thus large sepa- 
ration) of the graphene layers, this approximation of the bilayer potential as a superposition 
of single layer potentials is expected to be good, but, in any case, is certainly sufficient for 
the qualitative understanding presented here. 

As a basis for this Hamiltonian we take the combined eigenkets of the unrotated and 



rotated layers, i.e. { 0^^}, { <f>i% 2 ^}- One should note that since each SLG basis set by 
itself is complete on M 3 , this is generally an over-complete basis set. On the other hand for 
minimal basis methods, such as the p z tight binding method in which the basis consists of 
a p z atomic orbital centered at every site in the crystal, a bilayer basis set consisting of the 
combined eigenkets from each layer is clearly isomorphic to the usual basis set that would 
be employed. 
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B. Reciprocal space properties of the bilayer system 

Here we describe the reciprocal space lattices corresponding to the various real space 
lattices introduced in Section II. First a note of nomenclature; we denote the reciprocal lattice 
vectors corresponding to the unrotated (rotated) real space vectors &i and a 2 (Rai and Ra 2 ) 
by bi and b 2 (Rbi and Rb 2 ), while the reciprocal space lattice vectors corresponding to 
the the real space commensuration vectors ti and t 2 are denoted by gi and g 2 . We shall 
refer to this latter reciprocal space lattice as the bilayer reciprocal lattice. 

The vectors gi and g 2 are found from Eq. ffTB"]) to be 



3(3g 2 + p 



[(p + 3g)bi + 2pb 2 ] (20) 



3(3g 2 + p 2 j 

for the case where 5 = 1 and from Eq. (JUL to be 



[-2pbi - (p - 3q)b 2 ] (21) 



Si = ^^ [ - {p - q)hl + 2qh2] (22) 
& = 3^^[-29b!-(p + ? )b 2 ] (23) 

for the case 5 = 3. The Brillouin zones associated with each of these sets of primitive 
vectors, {b l5 b 2 }, {Rbx, Rb 2 }, and {gi, g 2 }, are shown in Fig. [3]for the twist bilayer (p, q) = 
(1,5). For convenience of exposition these Brillouin zones (BZ) will be referred to by the 
abbreviations UBZ (for the unrotated BZ), RBZ (for the rotated BZ), and BBZ for the BZ 
of the bilayer reciprocal lattice. 

These bilayer reciprocal lattice vectors determine a map by which k-vectors in the UBZ 
and RBZ are related to those of the BBZ (the usual so-called "folding back" condition of 
k-vectors). It should be emphasized at this point that there are three separate k- indices 
in the problem as it is formulated here. We have a k-vector in the BBZ which is a good 
quantum number for the bilayer Hamiltonian and eigenkets, but we also have the k-indices 
of the single layer basis used to solve the Hamiltonian at this k, labeled by kx and k 2 . To 
solve the Hamiltonian at k, the single layer basis then consists of all those eigenkets which 
map back from the UBZ and RBZ to the point k in the BBZ. 
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TABLE II: Structure of the mapping of special K points of the unrotated (U) and rotated (R) 
layer Brillouin zones to the special K points of the bilayer (B) Brillouin zone, for the case 5 = 1. 
The designation of the special K-points corresponds to that of Fig. [3J 
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-+K* B 


Kb 


^Kb 


mod (q, 3) 


= 2 


K v 


-^K* B 




^Kb 






K v 




K* v 


^K* B 






K R 


-^K* B 


K B 


^Kb 






«*■ 




Kb, 


^K B 



An interesting, and for the nature of the interlayer interaction crucial, relationship exists 
between the special K-points of all these BZ's: to each special K-point of the BBZ is mapped 
back one of the special K-point from the UBZ and one from the RBZ. The precise manner 
in which this happens depends in a rather complex way on the p, q parameters of the 
commensuration, detailed in Tables [TT] and IHIl The overall scheme, however, is clear: for 
commensurations characterized by 5 = 1 /^-points connected by AK (see Fig. |3]) map to 
the same special K-point of the BBZ while, for 5 = 3, it is the conjugate pairs (K v , Kr) 
and (Ku, K R ), separated by AK\ and AK2 respectively, that map back to the same special 
K-points of the BBZ (these special K-points and their separations are indicated also in 
Fig. E]). 

Thus without any layer interaction, i.e. what translational symmetry alone requires, 
is that the Dirac cones from the unrotated and rotated layers are mapped to the special 
K-points of the BBZ. With no layer interaction we therefore find two degenerate Dirac 
cones situated at each special K-point of the BBZ. It should be stressed that this mapping 
is particular to the K star; a similar map does not, for example, exist for the M star. 
Interestingly, this implies that M-point chiral fermions (which may be generated by the 
application of a periodic scalar potential^) would behave rather differently in this context. 

We may now consider what happens to this degeneracy when we turn on a layer inter- 
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TABLE III: Structure of the mapping of special K points of the unrotated (U) and rotated (R) 
layer Brillouin zones to the special K points of the bilayer (B) Brillouin zone, for the case 5 = 3. 
The designation of the special K-points corresponds to that of Fig. [3J 
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= 1 
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= 2 


mod (p, 3) = 1 


K v 


-+K* B 








K* v 


^K B 








Kjj 


Kg 










-+K* B 




-> K B 


mod (p, 3) = 2 











K fi K B K fi -> K B 

K* K B K* -+ K* 



action. In general, of course, such an interaction would result in a splitting of the Dirac 
cones, however this is not what happens for the case of mutually rotated graphene layers. 
The key to understanding this, as we now describe, lies in the remarkable behavior of the 
overlap elements of the bilayer Hamiltonian, Eq. ( 1T8|) . with states from the mutually rotated 
graphene layers. 



C. Matrix elements of the bilayer Hamiltonian 



Given a bilayer potential of the form V'- 1 -' + V^, and a basis set of single layer eigenkets 



{ ^iik;i)}'{ ^iak )}' ^ ne electronic structure will be determined by interlay er matrix ele- 



c(2) 



ments of the type 



iiki 



c(2) 



J i k /• Using this matrix element specific example, we 
now show how one may derive a general condition that determines whether such a matrix 
element vanishes or not. Using a plane wave expansion for each of the objects in this matrix 
element, i.e. 
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V« = J2 V$e< - (24) 
*2.M = £it + *G a Me i(k ' +RG *' (26) 



RG 2 



we find 



; iiki 



it) = E ( E / rf Ht +Gl+G ;(^ 

!,RG 2 J / 



Gi ,R,G 



E ^ki+G^ 2( 5ki+ Gl =k2+R G 2 ( 27 ) 
Gi,RG 2 



where we have made the convenient substitution G x — G 1 — G\. As the structure of Eq. ([2] 
arises simply from the differing in-plane translation groups of the constituent objects, any 
such interlayer matrix element may be cast into this form (although the coefficients Cki+G^ 2 
will obviously be different). 

Clearly, it is the Kronecker delta term in Eq. ( 1271) that is the most significant consequence 
of rotation. This Kronecker delta term ensures that in the double sum over Gi and RG2 
only terms satisfying 



d = RG 2 + k 2 - kx (28) 

contribute. This is just a commensuration condition between the reciprocal lattices of the un- 
rotated and rotated layers. However, in contrast to the real space commensuration condition, 
ai = Ra2, this involves not only the geometry via the R operator, but also a dependence 
on the single layer states through the term k 2 — ki. This removal of contributions from 
the interlayer matrix elements is a direct consequence of the mutual rotation of the layers, 
and can be seen as a destructive interference of the quantum states from each layer. It is 
now clear that the advantage of the approach deployed here is that we have separated the 
symmetry aspects of the problem, which generate a selection condition for the coefficients 
^ki+G^ 2 ' f rom details of the electronic structure which are contained in the actual values of 
these coefficients. 



1(3 



Continuity of wavef unctions and potentials in real space implies that these coefficients 
will decay to zero with increasing |Gi j2 | and will be largest for Gi j2 at or near the origin. In 
addition, the coincident points between the lattices Gi and RG2 + k 2 — k x become increas- 
ingly separated as 9 — > 0° (just as the size of the real space commensuration cell diverges in 
this limit, see Section II). Since these coincident points represent the only symmetry allowed 
contributions to the bilayer matrix elements, we see that for sufficiently small 9 there can 
be at most one contributing term, which occurs in the case of a coincident point close to the 
origin of reciprocal space. For example, if ki = k 2 in Eq. (|28|) then Gi = RG 2 = is a 
solution and, as the coefficient Cq is generally non-zero, so in turn will the matrix element 
will remain finite for all 9. On the other hand, in the case where all coincident points are 
sufficiently far from the origin, the matrix element will vanish for sufficiently small 9. (An 
illustration of these two cases is given in the left and right hand side panels of Fig. [5]). 

As we shall now show the term k 2 — kx in Eq. (128]) results simply in a shift from the origin 
of the commensuration lattice (i.e., the lattice of coincident points) that would be found for 
the case k 2 — ki = 0. The relation between the term k 2 — ki and this shift, which can 
be found by solving Eq. ( 1281) . thus plays a crucial role in determining which bilayer matrix 
elements will vanish. 

The solution follows by the casting of Eq. ( 128]) into a Diophantine problem, exactly 
as outlined in Section II for the real space case. This, and the solution of the resulting 
Diophantine problem, are described in Appendix B. Here we summarize the results with 
the solutions expressed in terms of the unrotated reciprocal lattice vectors as mibi +m 2 b 2 . 
One finds that for the case 5 = 1 two possible solutions given by 




and 




\+7T \ % I 



Q P \2h-l 



2 



While for the case 5 = 3 one finds 
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m = a-\ \ I - 7e I ; I (31) 

and 



m = a I + W 2 ' + f r 1 " 2 ' 2 ] (32) 



1 \ 2q J 7 U + g / 6 9 I 2Zi - Z 2 
Here a and /3 are arbitrary integers, 7 = gcd(3g + q, 3q — p), and l± and Z 2 are the integers 
that result when k 2 — ki is expressed in coordinates of the bilayer reciprocal lattice, i.e. 



k 2 - ki = Zigi + Z 2 g 2 . (33) 

(Note that since both k x and k 2 fold back, under translations by the bilayer reciprocal lattice 
vectors gi and g 2 , to the same k-point of the BBZ then their difference can be expressed 
as integer multiples of gi and g 2 . Hence in coordinates of the bilayer reciprocal lattice the 
difference k 2 — ki will always be integer.) Clearly, in all cases the solutions are of the form 

aGi + /3G 2 + AG, (34) 

with the important constant shift AG determined by k 2 — kx- 

It should be noted that these expressions provide only a partial solution to Eq. (1281) . The 
reason is for this is that while (mi, m 2 ) must be integer valued, the shift terms are obviously 
not integer valued unless, e.g., both l\ and / 2 are divisible by 6^/7 in Eq. ( f29l) . This absence 
of a complete solution is due to the fact that, as shown in Appendix B, Eq. (128]) results in an 
inhomogeneous simultaneous linear Diophantine problem (in contrast to the homogeneous 
problem of the real space commensuration) which is known to have no analytic solution. 
However, as we now demonstrate, this partial solution provides sufficient insight into the 
selection rule, Eq. (1281) . that several generic features of the bilayer electronic structure may 
be elucidated. 



D. Decoupling of the Dirac cones 

Here we shall prove that for all commensuration cells with Nq greater than some critical 
value the Dirac bands from the unrotated and rotated layers will be effectively degenerate 
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in energy. Thus there will be a fourfold degeneracy at the Dirac points of the bilayer 
band structure, with the Dirac bands themselves twofold degenerate. While this derivation 
demonstrates the existence of such a critical value, the numerical value of this parameter 
will depend on the coefficients Cki+G^ 2 ano - can ' °f course ; only be determined by actual 
calculation of the electronic structure. 

It should be stressed that this result requires only (i) mutually rotated and weakly in- 
teracting layers of hexagonal symmetry and (ii) the low energy spectrum to located at the 
vectors of the K star; it is, therefore, applicable to other contexts in which graphene may 
be created. 

Our approach is based on a perturbative treatment and the use of the selection rules for 
terms in the matrix element sums derived in the previous section. Since in the absence of 
any interlayer interaction we have two degenerate Dirac cones at the special i^-points of 
the BBZ (see Section III B), then the first order energy shift will be given by the secular 
equation of degenerate state perturbation theory. 

One should note that, without interaction, the degeneracy at and away from the Dirac 
points will be different: fourfold at the Dirac point and twofold away. In fact, as our 
considerations are entirely based on the in-plane translation groups, and not on point group 
symmetry arguments, then there will be no fundamental difference in how we treat these two 
cases, and for simplicity we consider here the case of a twofold degeneracy. The resulting 
secular equation is then 

(iBu SH a \ U\ = (1) (a\ 
\6H 2 i SH 22 J \a 2 J \a 2 J 
where the elements 8Hij are given by 



5R 12 
5H 21 
5H 22 



i> (1) 



V< 2 > 



; iiki 
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llkl I T'«2K2 
^(2) I ,(D 



4 (2) 
«2k2 



j(l) 

4 (2) 



V 
V 



4 (2) 
«2k2 

4 (D 
; iiki 



(36) 
(37) 
(38) 
(39) 



and with V = (V^ 1 ) + \^)/2. Here (z'i, ki) and (i 2 , k 2 ) are the k-vectors and band indices 
of the states from the Dirac cones of the UBZ and RBZ that map to the BBZ Dirac cone. 
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Using the approach of the previous section we can now determine when the matrix ele- 
ments involved in Eqs. fl36|) -f )37|) vanish. In particular, we know that for a sufficiently small 
misorientation 9 the vanishing or not of these matrix elements is governed solely by the shift 
term AG, i.e. by k 2 — ki. 

To determine AG we must therefore specify the difference k 2 — k 1; that is the difference 
between the k vectors of the two single layer Dirac cone states that map back to the same 
k- vector in the BBZ. If we consider the case 5 = 3 then, reading from Table [TTT| we find that 
the special .fT-point pair (K^, Kr), or its conjugate pair, is folded back to the same special 
-fT-point in the BBZ. Clearly k-point differences are unchanged by shifting all k-points by 
some 5k, and therefore k 2 — k x = — Kr. 

Expressing — K# in coordinates of the bilayer reciprocal lattice we then find — 
K R = (q- p)/(27)(gi + g 2 ), i.e. that h = l 2 = (q - p)/(2j) in Eq. flU). Using this and 
substituting into the shift term of Eq. (13T!) we find 



AG = 




We can easily ensure this is integer valued (as it must be) by the choice q = p(l + 2n), 
with n an integer. We thus conclude that both the commensuration reciprocal lattice vectors 
and the shift from the origin diverge as q — > oo, that is, as 9 — > 0°. The first order shift 
will therefore be negligible for all q that result in commensuration cells greater than some 
critical size, which will depend on the particular form of the C k ^t Gi 2 , i.e. on details of the 
electronic structure. 

If we reflect that as 9 -> then both AK X = \K* V -K R \ -»■ 2/3 and AK 2 = \K V -K* R \ -> 
2/3 while, at the same time, g = |gi| = |g 2 | — > (see Eq. 1421) . then it is not surprising 
that AG diverges. The behavior of AG in the case where 5 = 1 is, however, not so 
clear. In this case i^-points separated by AK map back to the BBZ .fT-points, and so 
k 2 — k x = K(7 — Kr; = — K.* R , and this, as may be seen from Fig. [3] goes to zero as 9 — > 0. 
In this case neither Eq. [29] or Eq. [30] yield an integer valued AG and so they may not be 
used. The Diophantine problem can, however, be solved numerically with the result that, 
as shown in Fig. [U | AG| diverges as 9 — > in this case also. 

The question then arises if higher order terms in perturbation theory may lead to a split- 
ting of the Dirac cones. In fact, it is easy to show that such terms may lead only to an equal 
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FIG. 4: (Color online) Shown is the shift term |AG| corresponding to the separations k2 — ki = 
- K R for p = 1, 2, 17 (5 = 3), and k 2 - ki = Ku - K R for p = 3 (6 = 1). Note that this shift 
diverges in all cases as 9 — > 0. 
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shift of both bands. This can be seen by an examination of the quantities involved in higher or- 
ders of perturbation theory. Let us first consider the calculation of the shift of the unrotated 
layer Dirac band. All terms in the perturbation expansion (which we do not need to consider 



explicitly) will involve matrix elements 



V 



« 2 k 2 



and 



i j Ci ), and the un- 



perturbed eigenvalues, which are just those of single layer graphene. Now, if we consider 
the shift of the rotated layer Dirac band we see that the relevant matrix elements are either 



the conjugate of those involved in the former case, 



i 2 k 2 



V 



,,(1) 



r 2 k 2 



l(2) 
«2k 2 



J iiki 



^iiki 
V( 2 ) 



V 



j(2) 
i 2 k 2 



; ilkl 



or 
Since 



are equal by the symmetry of the bilayer 
the unperturbed eigenvalue spectrum is again that of single layer graphene we immediately 
see that all terms in the perturbation expansion for the eigenvalue shift of the rotated and 
unrotated layers will be identical, and hence also the final energy shift. 

Thus it is only the first order term that can break the degeneracy of the Dirac cones from 
each layer and, as we have shown above, this is zero for all Nc greater than some critical 
value. Since it is known from ab-initio calculations^ that this degeneracy is already very 
small for the smallest cell Nc = 28, corresponding to (p, q) = (1,3), we can then conclude 
that for all commensuration cells the Dirac bands will be effectively degenerate, with exact 
degeneracy in the 9 — > 0° limit or incommensurate rotations. 

When higher order terms in the perturbation expansion are unimportant the enforced 
vanishing of the 1st order term therefore leads to a decoupling of the Dirac cones. Higher 
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FIG. 5: (Color online) Example of the relationship between k2 — ki and the shift of the lattice of 
solution vectors of the equation Gi = RG2 + ^2 — ki. Shown in the middle panel are those k-points 
that fold back to the T point of the commensuration Brillouin zone, for the case (p,q) = (11,31). 
The left and right panels display the reciprocal lattices Gi and RG2 + k2 — ki along with the 
solution vectors for the two different cases of k2 — ki indicated in the central panel. See Section 
HIE] for details. 




orders in perturbation theory, while unable to split the Dirac cones, may, as we shall see 
subsequently, lead to a suppression of the the Fermi velocity of the degenerate Dirac cones 
and non-linear band warping for very small misorientation angles. 

E. Rotation angle versus cell geometry dependence of the coupling of single layer 
states 

Given a relative rotation 9 between the two graphene layers there exists an infinite set 
of integer pairs (p,q) that, via Eq. (EJ), reproduce the rotation angle to arbitrary accuracy. 
This includes also incommensurate rotations, corresponding to the limit of diverging p and 
q. Each of the bilayer unit cells in this set have differing in-plane primitive vectors ti and t 2 
and reciprocal lattice vectors gi and g2. The symmetry allowed coupling of two single layer 
states (ii,ki) and (12,^2), governed by gi and g2, can, therefore, be dramatically different 
for very similar 6. As the electronic structure of the bilayer is determined by the coupling of 
single layer states through the interlayer interaction, this situation is counter-intuitive: one 
would expect electronic properties to depend smoothly on the rotation angle. 

Recalling the analysis of Section IIII Cj we note that the magnitude of the interlayer 
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coupling is determined by the selection rule, Eq. ( )28|) . for the coefficients of the Fourier 
expansion of interlayer matrix elements, see Eq. f )27|) . As has been discussed, the coefficients 
in such a Fourier expansion decay with increasing Gi j2 and are largest at or near the origin 
of reciprocal space. Now, the shift term of, e.g., Eq. (l3~Tj) describes a scale relation between 
the resulting shift of allowed G-vectors in the Fourier sum, and the k-vector difference of 
the coupled states k 2 — ki: states coupled by 2p/j(iigi + 2 2 g 2 ) correspond to a shift of 
— (zi,i 2 ). Thus it is not the gj that control how single layer states couple, but rather the 
vectors = 2p/'yg i . 

This is illustrated in Fig. Shown in the central panel are the UBZ, RBZ, and BBZ of 
the (p,q) = (11,31) twist bilayer (rotation angle 9 = 23.16°). Indicated are two different 
possible couplings of single layer states: (i) k 2 — k x = g 2 , and (ii) k 2 — k x = llgi = gj[ . 
The solution vectors of Gi = RG 2 + k 2 — ki corresponding to each of these cases are shown 
in the left and right hand side panels respectively. The single layer states coupled by g^ 
have contributing vectors close to the origin indicating a strong coupling, while the states 
coupled by g 2 have contributing vectors far from the origin, indicating a comparatively 
weaker coupling. 

(c) 

The extent to which states connected by g- dominate other symmetry allowed couplings 
depends on the nature of the decay of the Fourier coefficients C^i+g^ 2 ' i.e. on details of the 
electronic structure. For graphene, as we show in Section IV, the situation is of a coupling 
dominated entirely by the vectors k 2 — ki = nig^ + ra 2 g 2 c ^ with < < 1. Regardless 
of the particular details of the Fourier coefficients, allowing |tj| — > 00 and |gj| — > for a 
given rotation angle 9, beyond a certain point, introduces no new interlayer coupling to 
the bilayer system as all additional symmetry allowed couplings will have zero magnitude 
by the arguments above. Thus the ultimate smoothness of electronic properties with 9 
is guaranteed by the fact of the decay the Fourier coefficients Cki+G^ 2 ' a natural result. 

(c) 

Clearly, the faster the decay of these coefficients the greater the dominance of a few g) in 
the coupling of single layer states and, hence, the less the electronic structure depends on 
details of the real space cell. 

To determine a general form of the gjf , we first calculate g = |gj| which, for 5 = 1 is 
given by 
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' 7 sinfl/2 = -^Lcosfl/2. (41) 



and for 5 = 3 by 



3= ; 27 = -^sinfl/2 = -1 cos 0/2 (42) 

v/3(p 2 + 3g 2 ) >/3p 3g 

To determine = |g^| we then multiply (7 by the inverses of the pre-factors to the shift 
terms in Eqs. (I25j) - (j32]) . i.e. by 69/7, 6\/3p/7, 2p/7, and 6\/3<z/7 respectively (the factors 
of arise from the different (l\ — 2/2,2/1 — l2) T structure of Eqs. ( 130]) and (l3"2l) ). In this 
way we find that the relevant is given by 

9<c> = 1 7 !sin 5' (43) 

and that gj = g^gi with gj the unit vectors formed from the primitive vectors gj. As 
expected, these coupling vectors depend only on the misorientation angle. Interestingly, g^ 
may be expressed in terms of AK, i.e. in terms of the separation of pairs of special K-points 
(Ku,K R ), see Fig. as 

g (c) = b^_ AK . (44) 


In Ref. LJ it was noted that AK constitutes an energy scale of the twist layer physics. Here, 
in our more general lattice treatment, we see that this is also an important reciprocal space 
length scale, describing the coupling of single layer states by the interaction. 



F. Reduction in Fermi velocity of the Dirac cones 

Here we examine the impact upon the degenerate Dirac cones of the bilayer of the non- 
zero matrix elements. Let us consider the Dirac cone state from the unrotated layer situated 
at K[/ + 5k which, therefore, has eigenvalue 

e (1) = si|5k| (45) 

where we set fivp = 1, and s\ is a sign indicating whether the eigenvalue belongs to the 
electron or hole Dirac cone. From Section fill CI we know that this state will couple with the 
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two states from the rotated layer Dirac cone also having k-vector K.y + 5k. These states 
have energies 

e {2) = s 2 \5k- AK| (46) 

In this equation AK is a reciprocal space vector connecting nearest neighbor Dirac cones 
(for 9 < 30° that we consider here), see Fig. |3j Due to the translational symmetry, there 

(c) (c) 

may also be coupling with states at ~Ku + 5k — riig\ — n 2 g 2 which have energies 

eg* = s 2 \Sk - AK - n lg ? - n 2 g { 2 c) \ (47) 

fcl 

where n\ i2 are integers. Note the use the coupling vectors g,- , that determine which single 
layer states couple through the interlayer interaction, see Section IIII El Choosing S = 3 we 
find from Eq. HH that 



g[ c) = V3AK gl (48) 
gj° = V3AK g2 (49) 

with gi and g 2 the unit vectors of the bilayer reciprocal lattice. Using this one may then 
show that 

I AK + n lg i c) + n 2 g ( 2 c) \ = AKr} nin2 (50) 

where 

vl in2 = 1 + n i™2 + 3(ni - n 2 )(n 1 - n 2 - 1) (51) 

i.e. that the 9 dependence of |AK + riig^ + Ti^g^l is entirely through AK. We may then 
expand en}n 2 as 

(eSJ 2 = Sk 2 + AK 2 V 2 nin2 - 2SkAK Vnin2 cos0 niri2 . (52) 

In this expression = Z(5k, AK + n 2 g^ + ^g^), and 5k = |5k|. The overall eigenvalue 
shift may be written to second order as 
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Se = y { ^f— + ^f— \ (53) 

^ AX) + M AX) _ J?) y ' 

where a£*^ a an d a njn 2 are coupling constants between the different states (i.e., squares of 
overlap elements of the kind discussed in Sections III C and III D). Finally, from Eqs. ( I50IT53"]) 
may then be derived that this shift leads to a reduction in the Fermi velocity given by 

^ = ^ L) (1-^) (54) 

(SL) 

where vf is the Fermi velocity of the misoriented layers, v F the Fermi velocity of SLG, 
and a an overall coupling constant given by 

(y s l- _|_ (y Sl + 

a=J2 rain % niW2 . (55) 

ni n 2 ? 7nm 2 

Thus a > as required for Eq. [5H to actually describe a reduction in the Fermi velocity. 

The angle dependence contained in AK can, equivalently, be expressed via the moire 
periodicity (see Section II) D leading to the form 

v F = v ( F SL \l-(3D 2 ) (56) 
where B is a related coupling constant. 

n 

Note that while Eq. (1541) is of the same form as that derived by Santos et al. in Ref. 
the origin is somewhat different. Rather than use a continuum approximation, as deployed 
in Ref. we have retained the lattice physics which, in fact, only enters in the form of the 
'coupling' primitive vectors, Eqs. ( )48|) and (149]) . 



IV. TIGHT-BINDING ANALYSIS 



We now turn to tight-binding calculations of the twist bilayer structures elucidated in 
Section II. Given that the number of atoms in the real space commensuration cell diverges as 
the rotation angle 9 — > 0° the tight-binding method offers perhaps the only way of exploring 
this interesting limit; ab-initio calculations are certainly not practical. Here we shall employ 
the same tight-binding method deployed by Santos et al. in their continuum approach to 
the twist bilayer-, one of the so-called environment dependent tight-binding methods^. 
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FIG. 6: (Color online) Tight-binding calculation of the Dirac point splitting in (i) AB bilayer, 
indicated by light shaded points and (ii) 9 = 30° ±8.21 twist bilayers, indicated by dark shaded 
open squares (9 = 38.21°) and diamond symbols (9 = 21.79°). The ab-initio values for the AB 
and twist bilayer splitting are given by the horizontal dot-dashed lines. Full/dashed lines are 
calculations with the environment dependence of the hopping matrix elements switched off /on. 
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In the original article by Tang et. al the environment dependent parameterization for 
carbon that these authors proposed was checked against a database including, amongst other 
three dimensional lattices, diamond and graphite, as well as a one dimensional carbon chain. 
Unfortunately, graphene and graphene-based structures were not part of this dataset, and 
it is thus important to first verify the accuracy of this method for this case. A sensitive test 
of accuracy is provided by the Dirac point splitting of graphene bilayer structures, which 
may be quite large in the case of the AB stacked bilayer (the ab-initio value is 0.78 eV) 
and on the other hand rather small in the case of twist bilayers e.g. 7 meV for both the 
9 = 30° ± 8.21° twist bilayers^. This latter case entails a particularly sensitive test as, 
although both the 9 = 38.21° and 9 = 21.79° systems have the same Dirac point splitting, 
the crystal geometries are actually quite different^. 

In Fig. [6] are shown calculations of the Dirac point splitting for the AB bilayer, as well as 
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FIG. 7: (Color online) Tight-binding band structures for rotation angles of = 21.79° (1,3), 9.43° 
(1,7), 2.65° (1,25), 1.47° (1,45), displayed respectively clockwise from top left. Shown is the band 
structure generated by direct tight binding calculation (wide black lines) and that generated with 
SLG basis approach, indicated by light shaded (green) lines. For comparison in each panel the 
folded back band structure of single layer graphene is shown (dashed lines) . The numbering of the 
panels corresponds to that of Fig. [3l 




the two twist bilayers with 9 = 30° ±8.21°. Surprisingly, one finds that even the Dirac point 
splitting of the AB bilayer is not well reproduced; a much reduced interlayer separation 
is required to recover the ab-initio result of 0.78 eV. In addition, the splitting of the 9 = 
30° ±8.21° twist bilayers is also underestimated and, furthermore, is quite different between 
the 9 = 38.21° and 9 = 21.79° cases. 

Fortunately, this situation is significantly improved by switching off the environment 
dependence of the hopping integrals, in which case the method is simply the usual tight- 
binding scheme with distance dependent pairwise hopping matrix elements (see Ref. 22[). 
Given this, a reasonable agreement with ab-initio calculations may be found. For a somewhat 
reduced interlayer distance of 3.17 A(5% smaller than the nominal experimental interlayer 
distance of 3.34 A), we find ~ 7 meV for the twist bilayer splitting, which is in very good 
agreement with ab-initio data, and a splitting of 0.63 eV for the AB bilayer, which is less 
good but still reasonable. This interlayer distance is indicated by a dashed vertical line in 
Fig. [6j Our choice of calculation method is therefore the parameterization of Tang et ai, 
but with the environment dependence suppressed, and the interlayer distance set to 3.17 A. 
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While the tight-binding method extends the domain of direct band structure calculation 
far beyond that which may be achieved by the ab-initio approach, for the particular case 
of bilayer graphene (or more generally graphene stacks) one may do significantly better. 
In particular, if it is only the low energy band structure that is of interest, the single layer 
graphene basis introduced in Section III is much more appropriate that the full tight-binding 
basis set. In this basis the bilayer Hamiltonian is 



[H(k)]j lkli2k2 
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< ^ > «2k2^ / 



(57) 



where z'lki and are the band- and k-indices of states that fold back to k (a reciprocal 
lattice vector in the BBZ), and the superscript of the kets has the same meaning as in 
Section III, i.e., (l)/(2) refers to eigenkets of the unrotated/rotated layers. Eigenvalues at 
k may then be obtained by diagonalising the matrix consisting of all states iiki and 12^-2 
that fold back to k. 

Matrix elements in Eq. ( 1571) will involve both on-site terms and terms involving inter- 



layer hopping integrals 



V( 2 ) 



and 



dn>) 



V 



A 2 ) 



where V = ^(V^ 1 ) + V^ 2 ^). The matrix elements 



' k ) are equivalent to three center hopping integrals, 



and so may be set to zero, while the matrix element 



J iiki 



V 



t(2) 

12 k2 



may be evaluated as 



^ nvni 

where R ni and R n2 are vectors from layer 1 (unrotated) and layer 2 (rotated) respectively, 
the sum is over all atoms in the twist boundary primitive cell, a?* k is the ^-coefficient of 
the eigenvector corresponding to the i n k„ state from layer n, n z a directional cosine, t ppa 
and tppn distant dependent hopping integrals, and Nc the number of carbon atoms in the 
bilayer primitive cell. 

The advantage of this approach is that if it is only the low energy band structure that is 
of interest, then only low energy eigenkets of the unrotated and rotated layers are needed 
in constructing the SLG basis. In practice, the number of such low energy states required 
will depend on the strength of the interlayer interaction, which for graphene layers is weak 
and leads to a rather rapid convergence of the basis set, as shown in Fig. |HJ For the 
case of the 9 = 3.48° bilayer shown the maximum size of the basis set is 2048 states, and 
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FIG. 8: (Color online) Convergence of two nearly degenerate eigenvalues on the lower branch of the 
bilayer Dirac cone. The rotation angle of the bilayer is 6 = 3.48° (corresponding to (p, q) = (1, 19) 
- see section II). Eigenvalues are calculated at k = K + ^j(r — K) with K and T the k- vectors of 
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numerical convergence is reached at 24 states. A further advantage for the special case of 
twist bilayers lies in fact that from the selection rule analyzed in Section IIII CI one knows 
that a priori many of the matrix elements in Eq. (I57l) will be zero. Inspection of the 
numerical value of the interlayer matrix elements, Eq. f )58|) . shows that they are negligible 
for k 2 — = riigi +?^2g2 with rij > 1 (see Section fill El for a description of g^); a typical 
case is illustrated in the inset of Fig. El 

For actual calculations one may then utilize a truncated basis in which only a fraction 
of the actual matrix elements required need be calculated. This extends by more than an 
order of magnitude the number of carbon atoms that may be considered: within the SLG 
approach Nc = 26, 068 could be treated within the same time that, by direct tight-binding 
calculations, a system of Nq ~ 1200 could be calculated. 

We now consider the low energy electronic structure of a selection of twist bilayers, shown 
in Fig. [7J with the band structure plotted along the MKT high symmetry points path in the 
BBZ. In panels 1-4 are shown 4 twist bilayers in the set p — 1, q G odd Z with q — 3,7, 25, 45 
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FIG. 9: (Color online) Relation between misorientation angle of the bilayer and Fermi velocity 
damping; dashed (black) line in the main panel is the best fit to Eq. (|54p for data points with 
9 > 5°. 
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(misorientation angles of 9 = 21.79°, 9.43°, 2.65°, and 1.47° respectively). In panels 1-3 are 
shown band structures generated by both direct tight-binding calculation as well as the SLG 
basis outlined above; clearly these two approaches lead to identical results, as expected. A 
number of interesting features may be noted from these band structures. Firstly, as expected 
from the general analysis of Section III the Dirac point always decouples, i.e., there is no 
splitting of the degenerate Dirac bands from each layer, see Fig. [TUJ On the other hand, 
in agreement with ab-initio calculations^^, one observes that away from the Dirac point 
the bilayer band structure clearly shows a perturbation due to layer interaction, as may be 
observed in Fig. [7] by the strong band hybridization at the T point. This, as discussed in 
Section IIII C\ is due to the term k 2 — k x in the selection rule Gi = RG2 + k 2 — when 
k 2 — k x = the interlayer matrix elements never vanish leading to a strong coupling for 
these states. In addition, one notes a reduction in Fermi velocity of the Dirac cone as the 
rotation angle 9 — > 0°. 

In Fig. [9] is shown this Fermi velocity suppression as a function of the rotation angle of the 
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FIG. 10: (Color online) Splitting of bands at the Dirac point in misoriented graphene layers. Shown 
are tight-binding binding calculations for supercells generated with p = 1 and q an odd integer. 
First principles calculations for supercells generated by (p,q) pairs of (1,3), (1,2), (1,5), (2,3), (1,7) 
in order of increasing N are taken from Ref. 
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twist bilayer. Clearly, this effect is quite substantial, and for 9 = 1.47° (the smallest angle 
calculated) Dp is only 5% of the value of SLG. Interestingly, for 9 > 5° the tight-binding 
data is very well described by Eq. ( 154)) . in the small angle limit, however, the failure of 
any fitting of the form Eq. ( 154"]) indicates the importance of higher orders of perturbation 
theory. Clearly, the band structure in the 9 — > 0° limit is profoundly altered from that of 
SLG with, however, the degeneracy of the Dirac bands from each layer preserved, as was 
proved must be the case in Section III D. (One should note that, in contrast to the first 
principles calculation of Ref.—, we find a small reduction in Fermi velocity for the case of 
9 = 9.43°; such a difference presumably reflects the fact that parameterisation of Tang et 
alj22l. not explicitly designed for low energy graphene physics, has scope for improvement.) 

As has been mentioned, this preservation of the degeneracy is a striking illustration of 
the singular nature of the 9 — > 0° limit; for any finite 9 one has a fourfold degeneracy at 
the Dirac point, while at 9 = one has the AB bilayer with a twofold degeneracy at the 
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Dirac point, and with the other Dirac bands bands split by 0.78 eV. An interesting question, 
which we shall only pose here and not answer, is how this Fermi velocity reduction would 
be altered by both charge self-consistency and many-body effects. Both of these may be 
expected to become more important as the Fermi velocity is reduced, and may dramatically 
change the nature of the small angle electronic structure. 

Graphene stacks grown on the C-face of SiC typically have Fermi velocity reductions 
of 20-30% which then implies, assuming that such a reduction is entirely due to rotation, 
misorientations of 9 > 5°. Given that the formation energy of a twist bilayer increases as 
9 — > 0°—, with the minimum defect energy for 30° ±2.20°— , it makes sense that, on average, 
misorientation angles with 9 < 5° should play a less important role. On the other hand, 
it is clear that samples with 9 < 5° do exist; STM experiments detect moire patterns that 
correspond to angles in the 1.9° — 19° ranged In contrast, for a graphene slab dominated by 
30° ± 2.02° rotations, or simply large angle rotations (such a system was studied in Ref. |8|, 
see also Ref. [23[) one would expect to have, in addition to a Dirac spectrum over a wide 
energy range, a Fermi velocity exactly that of SLG. 



V. CONCLUSIONS 



To conclude we have given a complete description of the possible commensurations of 
graphene layers misoriented by some angle 9. We find that the condition for a commensu- 
ration to occur is that, expressed in lattice coordinates of the unrotated layer, the rotation 
matrix connecting the layers be rational valued, and thus the complete set of commensu- 
rations is described by two integers, which we denote {p,q)- For any such bilayer, we have 
shown that the K points of the unrotated and rotated layers map directly to K points of 
the bilayer Brillouin zone; a fact that plays an important role in the interlayer interaction. 

We have further shown that the nature of the interlayer interaction may be understood 
by a k-dependent interference condition, that may be expressed as a commensuration of the 
single layer reciprocal lattices. This guarantees the decoupling of the Dirac point and the 
degeneracy in the Dirac cones from each layer, but does not preclude interactions between 
all states. These latter interactions in fact give rise to a reduction of the Fermi velocity in 
the 9 — y 0° limit, and we find a form of this Fermi damping which agrees with that presented 
by Santos et a/.-, although our derivation is independent of any continuum approximation. 



33 



As an interesting consequence of this analysis, we are able to show that the bilayer electronic 
structure will, in general, depend only on the misorientation angle of the layers and not on 
the details of the real space unit cell. 

To complement this general analysis we have calculated band structures of a wide range 
of graphene twist bilayers via the tight-binding method. By the introduction of a basis of 
single layer graphene eigenkets, which we show to be significantly more efficient for the case 
of the twist bilayer, we are able to probe the band structure in the small angle limit with 
relative computational ease. We find Fermi velocity reduction that, for rotation angles in the 
range 5° < 9 < 30° agrees very well with the form presented here and in Ref. but that 
for 9 < 5°, where at 9 = 1.47° the Fermi velocity is only 5% of the SLG value, the reduction 
cannot be described in this way. In fact, the Dirac bands in the small angle limit, while 
guaranteed to be degenerate, show non-linear distortion away from the Dirac point. Thus 
the graphene twist bilayer encompasses a wide range of electronic behavior, from essentially 
SLG behavior for large angle rotations to quite different behavior in the small angle limit 
which, nevertheless, shares important features with the large angle case. While small angle 
rotations^ (as low as 1.9°) have been observed experimentally, the electronic properties of 
such low angle misoriented layers has yet to be experimentally explored. The possibility of 
a graphene type behavior different from both SLG and the AB bilayer makes interesting the 
further study of this low angle limit. 

The authors acknowledge Deutsche Forschungsgemeinschaft for financial support, and 
gratefully recognize collaboration within the Interdisciplinary Centre for Molecular Materials 
at the University of Erlangen-Niirnberg. 

Appendix A: Derivation of the commensuration lattice primitive vectors 

In this Appendix we wish to determine the primitive lattice vectors of the commensuration 
lattice given by Eq. (ITT]) , i.e., the lattice defined by 



We first notice that only coprime p, q correspond to unique solutions. Given this, a 
necessary condition for recovering lattice vectors is obviously the elimination of the greatest 
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common divisor (gcd) of the components of all the (integer valued) vectors in Eqs. (fTTj) and 
(TT2]) . This leads to following commensuration primitive vectors: 

tl .l(* + 3 'U = l( 2P | (A2) 
1 \ -2p J 1 \-p + 3q / 

with 7 = gcd(3g + p, 3q—p), and where we have used the fact that gcd(x, y) = gcd(x + cy, y) 
with c an arbitrary integer. Possible values of 7 involve the additional parameter 6 = 
3/gcd(j>, 3), and are displayed in Tabled! 

To prove that these are indeed the primitive vectors requires further that there is no 
linear combination of them yields integer valued vectors of smaller length. In fact, as we will 
now show, the vectors given in Eq. ( 1A2|) are primitive only for the case where 5=1, and 
that when 5 = 3 it is the linear combinations l/3( — t x + 2t 2 ) and l/3(— 2ti + t 2 ) that form 
the primitive vectors of the commensuration lattice. 

We first take a linear combination of the supposed primitive vectors as follows: 

|t 1 + |t 2 = t; (A3) 

where ai,Pi,N G Z and the index i = 1,2. By eliminating common factors of ojj and Pi 
from Ni we may choose gcd(o!j,/3j) = 1. Defining t = \t\\ = | "t 2 1 and t! = \t[\ = |t' 2 |, ti and 
t 2 are primitive only if there exists no c^, Pi, Ni such that t! < t. Suppose that for a given 
set of these parameters t' > t then we may write 

Nit < N^ = |oit a + Pit 2 \ < (N + \Pi\)t (A4) 

so that if t' > t then |aj| + > Ni and hence \ai\ + < Ni implies t' < t. Thus if 
there exist cti,f3i,Ni such that + < Ni and the linear combinations Eq. ( 1A3|) remain 
integer valued then the supposed primitive vectors ti and t 2 are in fact not primitive. 
From Eq. (|A2]) and Eq. (jA3]) we find 



\Oti 



2A 



iN 

(2oi + Pi 

iN 



-p 
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where Z\,z<i G Z in order that t- be integer valued. As gcd(di, fa) = 1 and gcd(p, q) = 1, 
then for z\$ to be integer valued as claimed we require jNi to be a common factor of the 
coefficients of p and q on the left hand sides of Eqs. (1A5I) and (]A6j) . Using the fact that 
gcd(«j, = 1, we find that possible values of gcd(3«j, aj + 2/3j) and gcd(3/3j, 2a« + are 
1,2,3,6. Therefore, as the minimum possible value of 7 = 1, then the maximum possible 
value of Aj is 6. As \cti\ + \fa\ < Aj this in turn restricts the possible values of a i} fa. 

We have used the requirement that t- be integer valued to place a condition 7Aj < 
6, not Ni < 6, and hence we need to eliminate (oti,fa,Ni) that lead to non-integer 
t[. From the full set of {(cKj, fa, Ni)} this then leaves only the cases (aii,fa,Ni) = 
{(1,1,3), (—1, 2, 3), (—2, 1,3)}, and using these latter two we finally find the new primitive 
vectors 

tl ^(-H,t 2 = i( 2 « ), (A7) 
1 \ 2q j "t\-q + qj 

However, for the case 5 = 1 then we have 3 | 7 and 3 | p and since gcd(p, q) = 1 then 3 { q 
and ti and t 2 given by Eq. ( 1A7|) cannot be integer valued. Hence for 5 = 1 then Eq. (IA2|) 
are already the primitive vectors of the commensuration lattice. On the other hand if 5 = 3 
then if both p and q are odd then 7 = 2 and ti and t 2 given by Eq. (1A7j) are integer valued 
while if one of p, q is even then 7 = 1 and again this is so. 

To summarize we find that for the case 5 = 1 the the vectors given by Eq. (|A2[) are 
already primitive, while for the case 5 = 3 instead Eq. ( ]A7j) gives the primitive vectors. 

Appendix B: Solution of equation Gi = RG 2 + (k 2 — ki) 

We wish to determine the solutions to the equation 

G 1 = RG 2 + (k 2 -k 1 ). (Bl) 

This represents a similar equation to the real space commensuration equation, but with an 
additional term (k 2 — ki). Utilizing the coordinate system of the unrotated reciprocal lattice 
we may write write Gi = mibi + m 2 b 2 and G 2 = nibi + n 2 b 2 where m = (mi,m 2 ) T and 
n = (ni,ri 2 ) r must be integer valued. Furthermore, the term (k 2 — ki) is integer valued 
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in the coordinate system of the bilayer reciprocal lattice, i.e., (k^ — ki) = Zigi + hg2- The 
transformation from the bilayer to unrotated reciprocal lattice coordinate systems is 



■BU 



(B2) 



-p + q —2q 
2q —p — q 

and the rotation operator, transformed to the unrotated reciprocal lattice coordinate system, 
is be found to be 



where, as before, we have 



(B3) 




H = 2pq (B4) 
i 2 = 3q 2 -p 2 (B5) 
2 3 = 3q 2 +p 2 . (B6) 
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Using these we may rewrite Eq. fIBlj) as 

7 I -P + Q - 2 ? 

2q —p — q 

Our solution of this equation is based on diagonalising R^ and Tbu which, it turns out, 
may be simultaneously diagonalised. We find the eigenvalues of R^ to be 

a ± 7-ET ( B =V 

p=F iVoq 

and those of Tbu to be 

The eigenvectors in both cases are given by 

u± = —\ 2y + V J (BIO) 



y/2 

Using these results we may rewrite Eq. flB7j) as 



1 
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U^m = (U^RzUQlT 1 !! + (U^Tcu^U-H (Bll) 
Equating the real and imaginary parts of this equation we find 

(n 2 + m 2 )p = ((2ni - n 2 ) - (2m 1 - m 2 ))q - 7Z2 (B12) 
((2m - n 2 ) + (2m! - m 2 ))p = (m 2 - n 2 )3q - 7 (2Zi - h) (B13) 

and by introducing the new variables n 3 = 2ri\ — n 2 and m 3 = 2m\ — m 2 we can recast these 
equations CIS db Diophantine problem: 

(m 2 + n 2 )p = (n 3 - m 3 )q - 7/2 (B14) 
(n 3 + m 3 )p = (m 2 - n 2 )3q - 7(2/1 - l 2 ) (B15) 

By absorbing the terms involving 7 either in the coefficient of p or q these equations may 
now be solved by inspection giving 



7/2 

m 2 + n 2 -\ 1 p = (n 3 - m 3 )q (B16) 

p 

p 



n 3 + m 3 + 1^—1 1 p = ( m2 - n 2 )3q (B17) 



which leads to the solutions 



B = I [" + « 2 « _JL S (B18 ) 



1 \ 2p J 7 \- p + qj 2p \ i 2 

and 



m=0 I + />i 25 -f "l (Bl») 



7 \ 2q I ^ \p + qj 2 P \h 
with a and /3 integers. Note that we have immediately written down the solution in terms of 
primitive vectors; which may be done in a similar fashion to the real space case (Appendix 
A). 

Alternatively we may write 
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(m 2 + n 2 )p = [ n 3 — m 3 — q (B20) 
(n 3 + m 3 )p = [m 2 -n 2 - — —J 3g (B21) 



which in turn leads to the solutions 




7 / /i — 2/ 2 
6 <? I 2h - h 



(B22) 



and 




(B23) 



7 / ^1 — 2/ 2 
6 P ^2/i - l 2/ 

The case where 5 = 1 proceeds in exactly the same manner, the only difference being a 
different form for the transformation matrix Tbu- The solutions are then found to be 

-^(T)^t + V-- I " ;i 

and 





l+f r 2/2 1 
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6 P I 2Zi - Z 
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6 P V 2Z X — Z 
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